Details about this script

On 05 April, 2024 I worked on(and if you are at the html knitted) this script.
The aim of the R Markdown script(or the output) is:
MPA Visualization - NS Germany
These files were downloaded from: https://sdi.eea.europa.eu/catalogue/srv/api/records/dae737fd-7ee1-4b0a-9eb7-1954eec00c65
Natura 2000 (vector) - version 2021 revision 1, Oct. 2022
🗺️ Zoom in to check the details in the OpenStreetMap in the background.
👆 Click on shapes to see German names.

General notes for the reader:

Please follow the table of content (outline) for keeping track of the steps. It must be on the left side.

DISCLAIMER:
Some code chunks could be eval=F to not to run the chunks but still include them in reporting, or results=‘hide’might hide some results to avoid cluttering.
Since the output document is shortened to demonstrate part of authors’ skill set.

Housekeeping

Check the needed packages:

##print Rversion:#####
R.version.string

## [1] "R version 4.2.2 Patched (2022-11-10 r83330)"

#important for rmarkdown####
knitr::opts_chunk$set(echo = T, eval = T, fig.keep="all", cache = T) #DEFAULT:echo = T, eval = T, #eval=F does not run the chunk
knitr::opts_knit$set(root.dir = "~//")
#getwd()


#housekeeping
##packages####
options(rlib_downstream_check = FALSE) #dont check package downstream
#require(sp)
require(sf)

## Lade nötiges Paket: sf

## Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE

require(rnaturalearth)

## Lade nötiges Paket: rnaturalearth

#library(knitr)
require(leaflet)

## Lade nötiges Paket: leaflet

#library(cowplot)
library(ggspatial)

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(knitr)
library(kableExtra)

## 
## Attache Paket: 'kableExtra'
## 
## Das folgende Objekt ist maskiert 'package:dplyr':
## 
##     group_rows

library(units)

## udunits database from /usr/share/xml/udunits/udunits2.xml

##print sessionInfo:#####
#sessionInfo() #for reporting on package versions #commented out for privacy

Paths, parameters…

LOADING

Loading files that were downloaded:

# Set the path to your GeoPackage file
gpkg_path <- paste0(shp_path,"EEA-Natura2000/eea_v_3035_100_k_natura2000_p_2021_v12_r01/GPKG/","Natura2000_end2021_rev1",".gpkg")

EXPLORE

# Read the GeoPackage file
st_layers(gpkg_path)

## Driver: GPKG 
## Available layers:
##            layer_name geometry_type features fields
## 1  NaturaSite_polygon                  27020      5
## 2           BIOREGION            NA    28829      3
## 3   DESIGNATIONSTATUS            NA    78751      5
## 4    DIRECTIVESPECIES            NA     2315      9
## 5            HABITATS            NA   152089     16
## 6        HABITATCLASS            NA   139836      4
## 7     NATURA2000SITES            NA    27027     23
## 8        OTHERSPECIES            NA   334086     13
## 9            METADATA            NA        4      2
## 10             IMPACT            NA   200695      7
## 11         MANAGEMENT            NA    32145     14
## 12            SPECIES            NA   388303     19
##                         crs_name
## 1  ETRS89-extended / LAEA Europe
## 2                           <NA>
## 3                           <NA>
## 4                           <NA>
## 5                           <NA>
## 6                           <NA>
## 7                           <NA>
## 8                           <NA>
## 9                           <NA>
## 10                          <NA>
## 11                          <NA>
## 12                          <NA>

n2kALL <- st_read(gpkg_path,layer = "NaturaSite_polygon")

## Reading layer `NaturaSite_polygon' from data source 
##   `/home/VTI_AD/oerey/03_Data/shapes/EEA-Natura2000/eea_v_3035_100_k_natura2000_p_2021_v12_r01/GPKG/Natura2000_end2021_rev1.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 27020 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 746551.8 ymin: 905053.4 xmax: 6506127 ymax: 5298655
## Projected CRS: ETRS89-extended / LAEA Europe

#README<- st_read(gpkg_path,layer = "HABITATCLASS")

st_crs(n2kALL)

## Coordinate Reference System:
##   User input: ETRS89-extended / LAEA Europe 
##   wkt:
## PROJCRS["ETRS89-extended / LAEA Europe",
##     BASEGEOGCRS["ETRS89",
##         ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
##             MEMBER["European Terrestrial Reference Frame 1989"],
##             MEMBER["European Terrestrial Reference Frame 1990"],
##             MEMBER["European Terrestrial Reference Frame 1991"],
##             MEMBER["European Terrestrial Reference Frame 1992"],
##             MEMBER["European Terrestrial Reference Frame 1993"],
##             MEMBER["European Terrestrial Reference Frame 1994"],
##             MEMBER["European Terrestrial Reference Frame 1996"],
##             MEMBER["European Terrestrial Reference Frame 1997"],
##             MEMBER["European Terrestrial Reference Frame 2000"],
##             MEMBER["European Terrestrial Reference Frame 2005"],
##             MEMBER["European Terrestrial Reference Frame 2014"],
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]],
##             ENSEMBLEACCURACY[0.1]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4258]],
##     CONVERSION["Europe Equal Area 2001",
##         METHOD["Lambert Azimuthal Equal Area",
##             ID["EPSG",9820]],
##         PARAMETER["Latitude of natural origin",52,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",10,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",4321000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",3210000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (Y)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (X)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Statistical analysis."],
##         AREA["Europe - European Union (EU) countries and candidates. Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Iceland; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal including Madeira and Azores; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain including Canary Islands; Sweden; Switzerland; Türkiye (Turkey); United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State."],
##         BBOX[24.6,-35.58,84.73,44.83]],
##     ID["EPSG",3035]]

#transform projection for leflet #4326 is the EPSG code for WGS84
n2kALL<-st_transform(n2kALL, crs = 4326) 

# Explore the structure of the spatial data
print(n2kALL)

## Simple feature collection with 27020 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -32.36806 ymin: 27.58768 xmax: 34.10149 ymax: 70.02586
## Geodetic CRS:  WGS 84
## First 10 features:
##     SITECODE                       SITENAME MS SITETYPE
## 1  PLH020044            Stawy Sobieszowskie PL        B
## 2  PLH020088                Dalkowskie Jary PL        B
## 3  PLH220055               Bunkier w Oliwie PL        B
## 4  PLH020051 Irysowy Zagon koło Gromadzynia PL        B
## 5  PLH020069                  Las Pilczycki PL        B
## 6  PLH020068         Muszkowicki Las Bukowy PL        B
## 7  PLH020070             Sztolnia w Młotach PL        B
## 8  PLH020078                  Kumaki Dobrej PL        B
## 9  PLH020084          Dolina Dolnej Baryczy PL        B
## 10 PLH020087         Gałuszki w Chocianowie PL        B
##                     INSPIRE_ID                           geom
## 1  PL.ZIPOP.1393.N2K.PLH020044 MULTIPOLYGON (((15.67494 50...
## 2  PL.ZIPOP.1393.N2K.PLH020088 MULTIPOLYGON (((15.87839 51...
## 3  PL.ZIPOP.1393.N2K.PLH220055 MULTIPOLYGON (((18.55123 54...
## 4  PL.ZIPOP.1393.N2K.PLH020051 MULTIPOLYGON (((16.35672 51...
## 5  PL.ZIPOP.1393.N2K.PLH020069 MULTIPOLYGON (((16.94962 51...
## 6  PL.ZIPOP.1393.N2K.PLH020068 MULTIPOLYGON (((16.95908 50...
## 7  PL.ZIPOP.1393.N2K.PLH020070 MULTIPOLYGON (((16.54941 50...
## 8  PL.ZIPOP.1393.N2K.PLH020078 MULTIPOLYGON (((17.21745 51...
## 9  PL.ZIPOP.1393.N2K.PLH020084 MULTIPOLYGON (((16.5485 51....
## 10 PL.ZIPOP.1393.N2K.PLH020087 MULTIPOLYGON (((15.71207 51...

# View basic summary information
summary(n2kALL)

##    SITECODE           SITENAME              MS              SITETYPE        
##  Length:27020       Length:27020       Length:27020       Length:27020      
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##   INSPIRE_ID                   geom      
##  Length:27020       MULTIPOLYGON :27020  
##  Class :character   epsg:4326    :    0  
##  Mode  :character   +proj=long...:    0

List of member states in n2kALL:

unique(n2kALL$MS)

##  [1] "PL" "EE" "MT" "SK" "DE" "PT" "DK" "SI" "RO" "SE" "IE" "HR" "FI" "CZ" "ES"
## [16] "FR" "IT" "LV" "AT" "LT" "CY" "LU" "NL" "BE" "BG" "GR" "HU"

#unique(n2kALL$SITETYPE)

filter out only Germany from n2kALL

#filterout only DE
n2kde <- n2kALL[n2kALL$MS == "DE", ]

#class(n2kde)

#plot what we have (takes time!)
#plot(n2kde)

rm(n2kALL)

cat("Total number of polygons after filtering:",length(st_geometry(n2kde)))

## Total number of polygons after filtering: 5200

SHAPE WRANGLING

make DE-NS bbox

# Define the bounding box based on your specified range
bbox <- st_bbox(c(xmin = mgmt_lon[1], xmax = mgmt_lon[2], ymin = mgmt_lat[1], ymax = mgmt_lat[2]), crs = 4326)

# Convert the bounding box to a simple features (sf) polygon
bbox_polygon <- st_as_sfc(bbox)

#plot(bbox_polygon)

leaflet() %>% 
  addTiles() %>%
  addPolygons(data= bbox_polygon)

get NS DE eez and territorial waters

# Set the path to your GeoPackage file
DEeez_path <- paste0(shp_path,"EEA-Natura2000/DEeez/")
DEeez_sf <- st_read(DEeez_path)

# Crop 'coast_sf' to the specified bounding box
DEeez_sf <- st_intersection(DEeez_sf, bbox_polygon)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

leaflet() %>% 
  addTiles() %>%
  addPolygons(data= DEeez_sf, color = "black") #%>%
  #addPolygons(data= DEeez_sf_buffered, color = "purple")

cut out bbox-DE-NS from n2kde

# Crop 'coast_sf' to the specified bounding box
n2kdeNS <- st_intersection(n2kde, bbox_polygon)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

leaflet() %>% 
  addTiles() %>%
  addPolygons(data= n2kdeNS, color = "black")

rm(n2kde)
print(paste("Total number of polygons after filtering:",length(st_geometry(n2kdeNS))))

## [1] "Total number of polygons after filtering: 208"

cut out EEZ shape from n2kdeNS

# Crop 'coast_sf' to the DEeez (but it also merges)
n2kdeEEZ <- st_intersection(n2kdeNS, DEeez_sf)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

leaflet() %>% 
  addTiles() %>%
  addPolygons(data= n2kdeEEZ, color = "pink")

rm(n2kdeNS)

print(paste("Total number of polygons after filtering:",length(st_geometry(n2kdeEEZ))))

## [1] "Total number of polygons after filtering: 43"

We want to focus on only the biggest of MPAs in german EEZ and territorial waters. I decided to calculate the area size and sort them accordingly.

sort shapes according to area size

#add column with area
n2kdeEEZ$area <- st_area(n2kdeEEZ) #Take care of units

class(n2kdeEEZ$area)

## [1] "units"

#check unit
attr(n2kdeEEZ$area, "units")

## $numerator
## [1] "m" "m"
## 
## $denominator
## character(0)
## 
## attr(,"class")
## [1] "symbolic_units"

#it is meter squares
#turn it into kmsquare by /1000000
n2kdeEEZ$area2 <- as.numeric(n2kdeEEZ$area/1000000)

summary(n2kdeEEZ$area2)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##    0.000    0.160    9.902  640.751  109.824 5289.394

n2kdeEEZ %>% 
ggplot(aes(x = area2)) +
  geom_histogram(binwidth = 10, fill = "black", color = NA) +
  labs(x = "Area", y = "Frequency", title = "Histogram of Area")+
  theme_minimal()

knitr::kable(as.data.frame(n2kdeEEZ) %>% 
               summarize(
                 TotalRows = n(),
                 RowsWithAreaGreaterThan100 = sum(area2 > 100),
                 RowsWithAreaGreaterThan50 = sum(area2 > 50),
                 RowsWithAreaGreaterThan40 = sum(area2 > 40),
                 RowsWithAreaGreaterThan20 = sum(area2 > 20),
                 RowsWithAreaGreaterThan10 = sum(area2 > 10)
               ))

TotalRows	RowsWithAreaGreaterThan100	RowsWithAreaGreaterThan50	RowsWithAreaGreaterThan40	RowsWithAreaGreaterThan20	RowsWithAreaGreaterThan10
43	11	15	15	17	21

PLOT

all polygons

leaflet() %>% 
  addTiles() %>%
  addPolygons(data= n2kdeEEZ, color="black")

smallestareas:<50 km2

smallestareas<- n2kdeEEZ %>% dplyr::filter(area2<50)
leaflet() %>% 
  addTiles() %>%
  addPolygons(data= smallestareas, color="black",
                  popup = ~SITENAME)

These must be partial shapes due to cutting out of the DE-EEZ and territorial waters.

biggestareas:>50 km2

n2kdeEEZ %>% nrow()
n2kdeEEZ %>% dplyr::filter(area2>50) %>% nrow()

biggestareas<- n2kdeEEZ %>% dplyr::filter(area2>50)

leaflet() %>% 
  addTiles() %>%
  addPolygons(data= biggestareas, color="black")

I decided to show only the shapes that are bigger than 50 kilometer square.

cat("Final number of polygons after filtering:",length(st_geometry(biggestareas)))

## Final number of polygons after filtering: 15

FINAL PLOT

check per site type

Since types overlap. Click to read the name. Use the legend on top right to choose site type.

unique(n2kdeEEZ$SITETYPE)

## [1] "B" "C" "A"

# split data
biggestareas.split <- split(biggestareas, biggestareas$SITETYPE)

# Define the color palette
pal <- colorFactor("Set1", unique(biggestareas$SITETYPE))

# Initialize Leaflet map
l <- leaflet() %>% addTiles()

# Iterate over each group in biggestareas.split
names(biggestareas.split) %>%
  purrr::walk(function(group) {
    l <<- l %>%
      addPolygons(data = biggestareas.split[[group]],
                  color = ~pal(biggestareas.split[[group]]$SITETYPE),
                  fillOpacity = 0.5,
                  fillColor = ~pal(biggestareas.split[[group]]$SITETYPE),
                  popup = ~SITENAME,
                  group = group)
  })

# Add layers control
l %>%
  addLayersControl(
    overlayGroups = names(biggestareas.split),
    options = layersControlOptions(collapsed = FALSE)
  )

save as ESRI Shapefile

n2KdeEEZbiggerthan50km2<-paste0(data_out,"/n2KdeEEZbiggerthan50km2")

# Write the sf object to a shapefile
st_write(biggestareas, n2KdeEEZbiggerthan50km2, driver = "ESRI Shapefile")

look at TABLE of final shapes

knitr::kable(as.data.frame(biggestareas) %>% dplyr::select(c("SITENAME","SITETYPE","SITECODE","area2")), 
             caption = "FINAL DE MPA SHAPES") %>%
  kable_styling(full_width = FALSE)

FINAL DE MPA SHAPES
SITENAME	SITETYPE	SITECODE	area2
Steingrund	B	DE1714391	173.72694
Borkum-Riffgrund	B	DE2104301	618.11929
Doggerbank	B	DE1003301	1676.77056
Nationalpark Niedersächsisches Wattenmeer	B	DE2306301	2574.97468
Sylter Außenriff	B	DE1209301	5289.39409
NTP S-H Wattenmeer und angrenzende Küstengebiete	B	DE0916391	4283.49474
Helgoland mit Helgoländer Felssockel	B	DE1813391	54.11029
Unterelbe	B	DE2018331	78.60430
Hamburgisches Wattenmeer	C	DE2016301	133.62144
Unterems und Außenems	B	DE2507331	51.94993
Schleswig-Holsteinisches Elbästuar und angrenzende Flächen	B	DE2323392	86.02672
Ramsar-Gebiet S-H Wattenmeer und angrenzende Küstengebiete	A	DE0916491	4305.92962
SPA Östliche Deutsche Bucht	A	DE1011401	3120.12815
Seevogelschutzgebiet Helgoland	A	DE1813491	1605.77646
Niedersächsisches Wattenmeer und angrenzendes Küstenmeer	A	DE2210401	3351.12014

This is just a plotting exercise. Please contact me for details: serra.oerey@uol.de 🧐

End of the document. - Serra Örey 💜

serraorey-MPAs-EU-DE

Serra Örey

Last compiled on 05 April, 2024